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Abstract 

A more accurate, stable, finite-difference time-domain (FDTD) algorithm is developed 
for simulating Maxwell's equations with isotropic or anisotropic dielectric materials. This 
algorithm is in many cases more accurate than previous algorithms (G. R. Werner et. 
al., 2007; A. F. Oskooi et. al., 2009), and it remedies a defect that causes instability 
with high dielectric contrast (usually for e ^ 10) with either isotropic or anisotropic 
dielectrics. Ultimately this algorithm has first-order error (in the grid cell size) when 
the dielectric boundaries are sharp, due to field discontinuities at the dielectric interface. 
Accurate treatment of the discontinuities, in the limit of infinite wavelength, leads to an 
asymmetric, unstable update (C. A. Bauer et. al., 2011), but the symmetrized version 
of the latter is stable and more accurate than other FDTD methods. The convergence 
of field values supports the hypothesis that global first-order error can be achieved by 
second-order error in bulk material with zero-order error on the surface. This latter point 
is extremely important for any applications measuring surface fields. 
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1. Introduction 

The Yee finite-difference time-domain (FDTD) algorithm |2] simulates electromag- 
netic waves in uniform, isotropic media with second-order error: i.e., the Yee algorithm 
simulates the frequency of a plane wave with an error that scales as 0(Aa:^) with grid 
cell size Ax. This paper presents a generalization of the Yee algorithm for non- uniform, 
anisotropic dielectric media, with particular attention to sharp transitions between dif- 
ferent dielectrics. (This generalization is also suitable for the intermediate cases, such as 
continuously-varying dielectrics, whether isotropic or not.) 

When different dielectric materials meet at a sharp interface, the discontinuity in the 
fields introduces greater error into the simulation. If the discontinuity is disregarded, 
operations such as field-interpolation typically have local 0(1) error at the interface (the 
error remains constant as Ax — because the field discontinuity remains constant). 
However, because the ratio of the cells cut by the interface to the total number of cells is 
0{Ax), the local error 0(1) is watered down by a factor of 0{Ax), leading to a "global 
error" of 0{Ax). (By global error, we refer to the error in a mode frequency or the 
average field error over an entire eigenmode, whereas local error is field error in a single 
cell. This relationship between global error and relatively high local error on a subset 
of cells has been proven rigorously for ID waves [3|, and demonstrated empirically for 
2D and 3D electromagnetic problems, e.g., with curved metal boundaries |4J and curved 
dielectric interfaces [SliHj.) 

When the dielectric constant varies continuously, the variation across a cell vanishes 
as Ax — 0, so it is not difficult to obtain local 0(Ax) error, hence global O(Ax^) error 

Of course we would like the same O(Ax^) error even with sharp dielectric boundaries; 
to the best of our knowledge, the first finite-difference approach to accomplish this was 
[5], which obtained global second-order error with local first-order error at the interface. 
Instead of considering the error of the discretized system for fixed frequency or wavelength 
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A and vanishing Ax, accuracy was demanded for fixed Ax and infinite A. Ref. fB| showed 
how to convert D to E exactly in the hmit of A — )■ oo. This led to global 0{Ax^) error. 

Unfortunately, the accurate method of [B] is unstable in the time-domain, because 
it uses an asymmetric inverse dielectric matrix (the linear operator that transforms the 
D field to the E field on the Yee mesh), which has complex eigenvalues, hence complex 
mode frequencies. While the imaginary parts of the frequencies are within the error 
[i.e., lm{u!) < 0(Aa;^/A^)], they lead to unphysical growth that becomes significant after 
sufficient time. Moreover, while well-resolved modes {Ax <^ A) may have slow growth, 
there are always modes with A ^ Ax which may grow quickly; thus machine-prccision- 
level noise eventually grows to overwhelm the desired signal. 

A symmetric inverse dielectric matrix was given by |5J, yielding an algorithm with 
0{Ax) global error (we will refer to this algorithm as "wc07"). This algorithm is stable 
at the low dielectric contrasts studied in ^ (where "dielectric contrast" is the ratio 
of dielectric constants between neighboring media in a simulation). However, we have 
recently found that wc07 is unstable for high dielectric contrast (see Sec. [5]), because the 
dielectric matrix is not always positive definite. 

In this paper we use the "triplets" concept of [6] to obtain stable algorithms in the 
time domain. If one can find symmetric and positive-definite (SPD) effective dielectric 
tensors acting on triplets of field components, then the tensors can be combined into an 
inverse dielectric matrix that is also SPD, which ensures stability in the time domain. We 
show that a small change to the wc07 algorithm [5] yields a new algorithm ( "wc07mod" ) 
with an SPD inverse dielectric matrix; wc07mod is therefore stable for arbitrary dielectric 
contrast. Using the same framework for stability, we provide yet another algorithm that 
is stable for arbitrary dielectric contrast; this algorithm simply uses the symmetrized 
matrix of [B]. The act of symmetrizing increases the error from 0{Ax^) to 0(Ax), 
but we find that this algorithm still has smaller error than the other effective-dielectric 
methods [3 Hill]. 

For example, for a dielectric contrast around 10, the new method reduces the fre- 
quency error by a factor of 2-3 at high resolutions, and it reduces the error of fields by 
a factor of 2-10. Reducing the error by a factor of 2 can reduce computation by a factor 
of 16, because the error scales as 0{Ax) while the computation time typically scales as 
0{l/Ax'^) (for a 3D problem). 

While examining the error, we show that for low dielectric contrast (less than 10 or 
so), the new algorithm yields very similar results to other effective dielectrics (OH] [7]), 
and the error is very nearly second-order up to high resolution. In other words, for low 
dielectric contrast, the error is dominated not by the field discontinuities, but by the bulk 
Yee algorithm. At sufficiently high resolution, however, we believe all these algorithms 
transition to first-order. As the dielectric contrast increases, the transition point moves 
to lower resolutions. 

We also examine the error in fields: the error at (or within a fixed number of cells 
of) the dielectric surface is 0(1). However, the field error is 0{Ax) at a fixed physical 
distance from the surface (N.B., that distance spans more cells as Ax diminishes). This 
supports the application of |3] to this problem; in other words, the local surface error is 
0(1), but only 0{Ax) cells are cut by the surface, so the global error is 0{Ax). 

The latter point may be very important for applications attempting to characterize 
surface fields with effective dielectric algorithms. Such attempts should be wary of errors 
in the surface fields, because the error does not decrease with cell size. However, the error 
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is probably small enough in many cases that it will not be a problem. And, the fields a 
fixed distance from the surface do become more accurate as the cell size is reduced. 

After a brief outline of algorithms discussed in this paper, we will present the dis- 
cretization of Maxwell's equations, reducing the problem of introducing dielectric to the 
problem of finding an inverse dielectric matrix S. Section [4] then describes how to create 
S to guarantee stability, assuming the ability to create local 3x3 effective dielectric 
tensors that are SPD (but otherwise unrestricted) , thus reducing the problem to finding 
the local effective dielectric. 

Sections [Sj [6| and [7] present three different methods for calculating local effective 
dielectric tensors. First we describe the effective dielectric that reproduces the wc07 
algorithm — this effective dielectric is not (always) positive definite, so it doesn't guarantee 
stability. Second, we modify wc07 slightly to guarantee stability. Third, we present an 
effective dielectric that guarantees stability, and has similar or better accuracy than the 
second option. All these methods require the same amount of computation for every 
time step. 

Subsequent sections present the error, in both frequency and field, of the different 
algorithms for different problems: 2D and 3D, isotropic and anisotropic, over a range of 
dielectric contrast from 5-100. 

2. Algorithms in this paper 

The following algorithms will be discussed in this paper: 

• wc07: the algorithm recommended by fSj [therewithin called variant (c) + (e)] — it is 
unstable for high dielectric contrast; 

(Rcf. [7, used wc07, with the improved dielectric averaging of [T] instead of [5] . 
These two averaging methods are identical for isotropic dielectrics; accuracy for 
anisotropic dielectrics may improve, but the error still converges as 0(Aa;). When 
we refer to wc07, we mean with the averaging of [1], but everything we say would 
apply just as well to [8].) 

• wc07mod: a stable algorithm, almost the same as wc07; 

• the "new" method: a stable and more accurate algorithm, but still with 0(Ax) 
error (sections |4] and [t]) ; 

• the second-order method of [6]: the only finite-difference algorithm with second- 
order error is unfortunately asymmetric, rendering it unstable for time-domain use 
(but still good for frequency-domain eigensolvers). 

3. The basic algorithm 

We want to simulate the dynamic Maxwell equations with dielectric: 

f = -VxE (1) 

f - ,2, 

E = iix,y,z)T> (3) 
4 



Figure 1: (color online) Field components in one grid coll of the Yee mesh. If the above cell has the 
3-integer index {i,j,k), then the labeled field components have the same index, e.g., Eijf^^. The D 
components are collocated with the corresponding E components. 



where £, = e ^ is the inverse dielectric tensor (which can vary in space). 

To discretize Maxwell's equations for computational work, we follow Yee To 
analyze the discretization, we treat the fields as large vectors, with each component 
representing a field at one point of the Yee mesh, shown in Fig. [l] (Discretization in 
time is irrelevant to this paper, so we retain continuous time derivatives.) 

For convenience, we label each component of a field vector with 4 sub-indices: e.g., 
the component Eijkx represents the electric field in the x direction at its Yee location 
in cell {i,j,k). The differential operators become matrices, with rows and columns each 
indexed by 4 sub-indices: e.g., one element of a matrix M is M(^ijk^-) (^iijikt^y 

We will not review the Yee discretization, since |5] details the relevant aspects, instead 
taking it for granted that the matrices C and represent the curl operators (the matrix 
representation of the curl of E is the transpose of the curl of B). We depart from Yee 
when we introduce the matrix S to discretize the linear relationship between E and Z?, 
yielding: 

d 

-Q-^B^jUt, = -[C£^]ijfc^ (4) 
d 

i^AjTcM = [C B]ijkti (5) 

Eijkfj, [^B]iji^f^ ^ ^ ^{ijk^),{i' j'k' y)E^i' j'k'v • (6) 

i' j' k' V 

These equations are a prescription for advancing the fields in time: B is advanced by a 
short time (using E)^ then D is advanced (using B), then E = ED, and the cycle repeats 
(in practice, this can be implemented with only two fields, E and B). Combining these 
into one equation 

Q^Bijk^ = -[CSC B]ijk^ (7) 



shows that the eigenvalues uP' (frequencies squared) of the —d^/dt^ operator are the 
eigenvalues of the CEC^ matrix. 

Simulating dielectrics therefore reduces to the determination of S such that 

1. S accurately represents the inverse dielectric tensor, ^{x,y,z), and 

2. CEC^ is diagonalizable with only real, non-negative eigenvalues. 

The first point addresses accuracy. The second addresses robustness/stability: if CSC'^ 
has negative or complex eigenvalues, some frequencies lu will be complex, and some 
fields will grow exponentially (unphysically) , ultimately overwhelming the simulation 
with noise. 

If S is symmetric and positive definite (SPD), then the second point above will be 
guaranteed: all modes will oscillate (with real frequency) without growing (c.f., [S], or a 
standard linear algebra text such as [5]) — at least for a sufficiently small time step. (As 
stated, temporal discretization is outside the focus of this paper). 

This paper is devoted to finding an SPD matrix !B that accurately represents dielectric 
media. 

4. Creating a stable H matrix from local ^ tensors 

Our approach to dielectric simulation falls under the "effective dielectric" category, 
specifying the matrix S to find E = 'E.D, while advancing D (and B) in time according 
to the unaltered Yee algorithm [5]. 

As in Ref. [6], we demand that S have the following properties. 

1. In case of uniform (anisotropic) dielectric, S involves only centered interpolations 
of fields (i.e., interpolations with 0{Ax'^) error for continuous fields). 

2. Within uniform, isotropic dielectric, S reduces to a multiple of the identity. 

3. S must be symmetric. 

4. S must be positive definite. 

The first property yields a local 0{Ax^) error within uniform dielectric (or even continuously- 
varying dielectric [5]). The second is for convenience and minimalism — we can still use 
the plain Yee algorithm for fields within any bulk isotropic region. Together, the third 
and fourth properties guarantee stability. (For comparison: wc07 satisfies 1, 2, and 3.) 

Before defining the S matrix, we need to introduce notation for an effective dielec- 
tric tensor involving field components of a single Yee cell. For example, we consider 
the "triplet" of three i?^ values labeled in Fig. IT] and the three I?^ values at the same 
locations, along edges that touch a common node(corner) of the cell. If the node index is 
{ijk), then this triplet comprises {Eijkx, Eijky, Eijkz), and similarly {Dijkx, Dijky, D^jkz)- 
(We call these components a triplet, instead of a 3-vector, because they are not collo- 
cated.) The triplets for E and D can be related by a 3 x 3 tensor, Ci^^^: 

Ejjfe = ^tjk^^ijk or S.jfc^ = ^ Ctjk^i^I^ijku (8) 

(the meaning of the + + + superscript will be explained shortly). Defined thus, S,tjk^ 
an effective (inverse) dielectric tensor for this particular triplet of E and D values. 
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Touching the same cell node are seven other (eight, in all) geometrically-identical 
triplets. Above, we chose a triplet with edges extending positively in each direction from 
their common node, but (due to the symmetries of the Yee mesh) we could equally well 
have chosen, e.g., the a;-edge extending in the negative x direction, namely (-B(i_i)jjti:, Eijky, Eijkz) 
and {D(^i_i)j).x,Dijky, Dijkz)- With this triplet, we associate a different effective dielec- 
tric, S^j^'^ ^ with the minus sign indicating that the and Dr^ edges extend in the 
negative x direction from node (ijk). 

We will construct (the large matrix) S from the (small) S.^^^ tensors, through an 
intermediate stage involving (large, block-diagonal matrices) 2='=='=='=. Symmetry and pos- 
itive definiteness transfer easily from one stage to the next: we will show that S is SPD if 
all the 'Ci^fc ^ Eire SPD; since the latter are mere 3x3 matrices, evaluating their positive 
definiteness is easy, numerically if not analytically. 

We start by creating block-diagonal matrices S+++, S S , where the 

blocks are the ^fj^^] for example, each block of S is S.^^^^ for some node {ijk). 
Precisely we define 

"+++ = A, 

"(ijkti)(i'j'k'y) — "Wk){i'j'k')<iijk,fii^ 

^(ijkn)(i'j'k'u) = ^i3,k+S^^,^lu^(ij.k+S^^),{i'j',k'+S^^) 

1- = ^ X 

"(iikii)(i'j'k'u) — 'ii.,j+S^y,k,tii^''{ij+Sf,yk),(i'j'+5^y,k') 

■^H = X (9) 

^{ijkfj.)(i'j'k'u) — 'ii.j+S^y,k+S^^,tj.u"(ij+S^y,k+5^^),{i'j'+S^y,k'+S^^) 

^{ijkp,){i'j'k'i') = ^i+S^,^,j+Sf,y.k+S^^,^lu^(■i+S^^:j+S^y,k+Sf,^),{i'+S^^J'+5^y,k'+S^^) 

where Si^i^j^'^^iijiy) is the Kronecker delta, equal to one when i = i', j — j', and k = k' , 
and otherwise equal to zero. A block-diagonal matrix is SPD if and only if each block is 
SPD; therefore, if each ^+++ is SPD, then S+++ is SPD, and likewise for all 
We finish by taking the S matrix to be the average, 

E^l + + + + E+- + E-+- + E-+ + E ) . (10) 

8 

Since the sum of SPD matrices is again SPD, S is SPD. Figure [2] illustrates how this 
scheme determines Eij^x from neighboring D values. 

Thus we can create a stable algorithm independent of the details of the ^ , which 
can be chosen to improve accuracy at the dielectric interface, within broad constraints: 
the ^^^^ must be SPD. In addition, when node (ijk) is within uniform dielectric and far 
from a dielectric interface, Cj^^^ must equal ^{x,y,z), where {x,y,z) is the position of 
node (i,j,k). This guarantees that, within isotropic dielectric, each S±±± is a multiple 
of the identity, and so S satisfies requirement 2, above. Furthermore, within uniform 
dielectric, the algorithm becomes identical to wc07 (see Fig. l3|, which uses centered 
interpolation to yield 0{Ax^) error (within uniform dielectric [5]). Therefore, S satisfies 
requirements 1 and 2 (as well as 3 and 4). 

It is easy to forget, when viewing the algorithm as a way to find a single E^ (as in 
Fig. g, that C+fcX' ^tjk^y^ ^"^^ ^tjk^z must come from the same (SPD) tensor £,±-+ . 
Indeed, it was forgetting which £,ijk,^v had to be mathematically related to each other 
that led to the instability of wc07. 
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Figure 2: (color online) A pictorial representation E = SD, showing how Ex in cell k) is found from 
neighboring components of D. For each of 8 "triplets" a local effective inverse dielectric ^ (a 3 X 3 matrix) 
converts D to E. Ultimately E^ is found using 8 different ^ matrices, one for each triplet involving E^- 
By averaging over all 8 triplets, Ex depends symmetrically on its neighboring Dy and D^, which yields a 
centered algorithm that, in uniform (or continuously-varying) dielectric, has second order error [5|. We 
note that 4 of the triplets use ^-^d 4 use j.- 
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Figure 3: (color online) A pictorial representation oi E = 'ED as suggested in [5]. in cell (i,j,k) 
is found from at the same location and the nearest Dy and values. The are interpolation 
operators (see [5]). 



To reiterate: each 5='=='=='^ matrix is block-diagonal, with 3x3 blocks Hfj^^), each of 
which represents the effective dielectric tensor around node {ijk). As long as the Cj"*"*^ 



are SPD, the 5='=='=='= are SPD. This further implies that the average, Eq. 10 is SPD. 

We have thus shown how to find a stable matrix S, given the ability to find 3x3 effec- 
tive inverse dielectric tensors that map a triplet of neighboring components {D^jDy, D^) 
to the {Ex, Ey^Ez) at the same locations. Moreover, within uniform (isotropic or anisotropic) 
dielectric, this algorithm is identical to wc07 (hence it satisfies requirements 1 and 2). 

It remains to find the effective dielectric tensors CijJ* that will accurately represent 
the real dielectric. Reference [HI showed how to find with local 0(Aa;) error, 

yielding global 0(Aa;^) error; unfortunately, those ■Ci^fc * are asymmetric. 

In the next section, we will describe the effective dielectric tensors that yield the 
wc07 algorithm; some of those tensors may not be positive definite, so that algorithm 
can be unstable. We then describe a modification to make it stable. However, another 
effective dielectric turns out to yield lower error (Sec. [7|. We have tried several other 
effective dielectrics, and found them less accurate, but always yielding ultimate 0(Ax) 
global error (even stairstepped dielectrics have 0{/S.x) error). 



5. The wc07 method, unstable at high contrast 

Experiment tells us that wcOT (the algorithm of [5]) can be unstable; therefore, we 
need concern ourselves no more with that algorithm. However, it is interesting to see 
exactly how these algorithms differ, so in this section we describe wc07 within the frame- 
work of this paper, and show why it does not meet the previously-described (sufficient, 
but not necessary) conditions for a stable algorithm. 

Figure |3] depicts the wc07 method recommended in ^ [where it is called (c)-l-(e)]. To 
find Ex from its nearest-neighbor values: 

1. perform a centered interpolation of Dy and Dz to the nearest cell node; 

2. find Ex[x) — £,xxDx, using the effective dielectric ^(^'^^ at the x-edge center; 

3. find Ex(y) — s}xy Dy at each node, using the effective dielectric at that node; 
similarly, find Ex(z) = s!xz Dz at each node; 

4. interpolate the two Ex[y) from each node to the center of the a;-edge where Ex and 
Dx are located; interpolate Ex{z) similarly. 
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Finally, is the sum of parts coming from E^(^^-^, ^x(y)^ E^^^y 
In the framework of the previous section, wc07 is a choice of effective 

^ijk.xx '^(i+l)jk,xx ^ijk.xx V / 

^ijk,vy ^i(j + l)k,yy 'iijk,yy 

C±+± _ f±±- — <r(^'=) 

^ijk^zz ^ij(k+l)^yy ^ijk^zz 

The tensors djfe^ ^Ijfe^ ^"^^ ^ijfe^ ^""^ ^'^ found from the averaging method 

of [T], where is the "average" of ^{x,y,z) over a cell volume centered at the node 

(lowest corner) of cell (i, j, k), and Clj^'' is the "average" over a cell volume centered at 
the location of Eijkx (the x-edge-center) , etc. 



These Cij^ ^ do not satisfy the conditions for the effective dielectric required by the 
previous section. E.g., S.^j'^^ is symmetric, but it is not necessarily positive definite. 
The reason is that the diagonal and off-diagonal elements of Ci^fe^^ come from different 
tensors: ^j"2> ^i^k^ ^ ^^Ak ' ^^'^ ^iik ^ ■ Each of these four tensors is SPD (using the averaging 



method of [T], c.f. Appendix B), but there's no guarantee that a tensor with a mixture 
of elements from those tensors is positive definite. 

Indeed, we have found experimentally that the wc07 algorithm can yield an instability 
when the dielectric contrast is high enough. We hasten to point out that we have used 
that algorithm successfully on a wide range of problems without noticing any instability. 
An instability seems to be more likely for higher contrast, and for larger and more 
complicated dielectric shapes. 

By using the Gershgorin circle theorem to place a lower bound on the eigenvalues of 
the S matrix (if the lower eigenvalue bound is positive, then S is positive definite, assum- 
ing S is symmetric), we can prove for many particular simulations that the algorithm is 
in fact stable. For example, we usually find that simulations with e < 10 are provably 
stable (on an individual basis, by examining S with the Gershgorin circle theorem). 

Many simulations appear stable for long times even when not provably stable. Of 
course, it's hard to know whether there might be unstably-growing modes that would 
dominate the simulation if run 100 times longer; and such unstable modes might interfere 
with precision measurements well before they become obviously apparent. 

For a 2D simulation of photonic crystal modes of a square lattice of isotropic dielectric 
discs in vacuum (radius r = 0.37a, where a is the lattice constant), we have seen that at 
e = 60 we can run simulations with 48^ cells for a time 3000a/c (where c is the speed 
of light) without seeing an instability (which rules out instabilities growing faster than 
7 ~ 0.01c/ a, which corresponds to growth of 16 orders of magnitude in 3000a/c). For 
64^ cells, however, an instability grows as exp(7t), where 7 « 0.5c/a. 

For the same problem, with contrast e = 100, it's harder to find stability for any 
resolution. For 32^ cells, an instability grows with 7 w 3c/a, and for 64^ cells, 7 « 6c/a. 



10 



6. wcOTmod: a small change yields stability 



In this section we make a small change to the wc07 algorithm that renders it stable. 
Although this algorithm, which we will call "wc07mod," is not the most accurate, we 
present it because it is a relatively trivial modification of wc07, and the resulting degra- 
dation of accuracy is interesting, in light of the small modification, which still uses the 
same averaging method to find the effective dielectric within a given cell-sized volume. 

In the language of this paper, this algorithm is described simply as 

^±±± ^ t±±± ^ ^(") ^2) 

(n) 

where £/^l is the "average" inverse dielectric tensor for a cell volume centered at the node 
of cell (i, j, k) — where averaging is done according to [1]. 



Only the diagonal elements of the effective dielectric change, compared to Eq. (11). 
In the language of we need simply replace, e.g., in Eq. (27e) of 5j or in step 2 of 
the wc07 algorithm: 

1 



and similarly for the yy and zz elements. 



f (") I 
^ijkxx ' ^{i-\-l)jkxx 



(13) 



7. The new method 

The most accurate local effective dielectric, from [B], yields local 0{Ax) error, but is 
unfortunately asymmetric (except for a few surface cuts: e.g., when a planar surface is 
parallel to a grid plane). Simply symmetrizing it increases the error to 0(1), but turns 
out to be more accurate than other (symmetric) effective dielectrics. 

Section |4] reduced the problem of finding a stable S matrix to the problem of find- 
ing SPD 3x3 matrices, e.g., C-jjT^, that map a triplet of neighboring components 
(I?j;, Dy, Dz) to the {E^, Ey, E^) at the same locations — in a way that accurately repre- 
sents the real dielectric. In this section, we describe the best such recipe that we have 
found. 

While we focus on finding the effective dielectric for a single triplet, we will omit the 
(H h) and (ijk) super- and sub-scripts, which identify the triplet. 

This effective dielectric, which is not a volume-averaged dielectric as in [51 [Tj, derives 
from Ref. [B], which finds the unique 3x3 tensor ^acc that guarantees that {E^, Ey, EzY" — 
£,iicc{Dx, Dy, -Dz)"^ will be exactly accurate in the limit of infinite wavelength and planar 
interface. In other words, ^acc will convert (D^, Dy, Dz) to {Ex, Ey, E^) with no error, 
given that the triplets are from the finite-difference (or rather, finite integration) repre- 
sentation of an infinite- wavelength solution of Maxwell's equations. Unfortunately, as we 
have mentioned, ^acc is not symmetric. 

To achieve stability in the time-domain, we will use 

?cff = \ (Cacc + Cc) ■ (14) 

This will be stable; its accuracy will be evaluated empirically. 

Finding ^acc is a lengthy process fully described in [B] , so we present only a terse recipe 
for converting a triplet {D^, Dy, Dz) to {E^, Ey, Ez) in the presence of two dielectric 
regions, ei and £2 (both symmetric tensors). 
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1. Within a small region around the triplet (we use the cell volume centered at the 
nearest node), the dielectric interface is nearly planar; find the unit surface normal 
n. 

2. Each electric field component E^^ is associated with a cell edge L^; for each edge, 
determine the fraction of its length in dielectric ei . See [6] for explicit definition 
of Lfj, (and in the following). 

3. Each component is associated with a dual-face area (centered at the Yee 
location of D^, perpendicular to /i); for each area, determine the fraction of the 
area in ei. 



Form the 3x3 matrices (for p — 1,2) 
^ [nn^](/- 
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(15) 
(16) 

(17) 
(18) 



where [nn-^] is the dyadic matrix with elements [nn-^J^i^ = fififiij, and / is the 
identity. 

5. The accurate effective (inverse) dielectric tensor is: 



We then symmetrize that to find 



: (Cacc + face) 



(19) 



(20) 



6. 



The above can fail if 11 is not invertible, and if the resulting ^cff is not positive defi- 
nite. Failure is ruled out for isotropic dielectrics [5], and for anisotropic dielectrics, 
it has not yet happened in our experience. Nevertheless, it's important to guard 
against pathological cases. We suggest checking every foff for these two problems; 
if one should occur, then substitute the effective dielectric from [T] (which is proven 
suitable in Appendix B). This will happen so rarely, if ever, that the global error 



will not be significantly affected. 



8. Simulation Results 

We tested various FDTD algorithms on different dielectric problems: a square lattice 
of 2D isotropic and anisotropic dielectric discs in vacuum; a square lattice of 2D vacuum 
discs (holes) in isotropic dielectric; and a cubic lattice of 3D dielectric spheres in vacuum, 
for both isotropic and anisotropic dielectric. For many of these cases, we also tested 
different dielectric contrasts. We define a to be the lattice constant, and Na the number 
of (square or cubic) cells per lattice constant, hence Aa; = a/Na- 
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Ultimately, the FDTD algorithms all show first-order error in frequency; the error in 
a mode frequency falls as 0{Ax), or 0{Ax/X), where A is a characteristic wavelength of 
the mode, with decreasing cell size Ax. However, at coarse resolutions (large Ax), the 
error often falls as O(Ax^) for low dielectric contrast. This may explain why previous 
studies have concluded that methods such as wc07 have second-order error — they did not 
explore high-enough resolution or contrast (of course, in practice, one may often reach a 
tolerable error level within the second-order regime, in which case the ultimate order of 
error may be irrelevant). 

The error convergence in surface fields was the same as in frequencies when we con- 
sidered the surface fields a fixed distance (e.g., a/8) away from the dielectric boundary. 
However, the error in fields a fixed number of cells (e.g., 3 Ax) away from the boundary, 
is 0(1), not 0{Ax). This supports our assertion that the local error at the boundary, 
due to discontinuous fields, is 0(1), but the global error is 0{Ax) because the ratio of 
boundary cells to total cells is 0{Ax). 

We performed the FDTD simulations with Vorpal [4 using the EDM method (TD] 
to extract accurate mode fields and frequencies. We compared these results with the 
frequency-domain algorithm of JiJ, which was shown to have second-order global error. 
For 2D simulation frequencies, we extrapolated results from Na = 512 and Na ~ 1024 
assuming second-order convergence to get a normative value with approximately 0{Ax^) 
error. 

We will show the most detailed convergence results for the "new" algorithm (Sec. [t]) 
for 2D anisotropic discs. Isotropic and 3D dielectrics show similar convergence, so we 
present only a few examples. 

We will show that the other FDTD algorithms, wc07 and wc07mod generally have 
similar or greater error compared to the new algorithm; and the other examples (3D and 
isotropic) show qualitatively similar convergence. 

For comparison, we also show frequency convergence for the second-order (but un- 



stable in the time-domain) method of "B] in Appendix A 



8.1. Convergence: 2D anisotropic discs 

We simulated TE modes (with — 0, = = By, and no variation in z) in a 

square lattice (lattice constant a) of dielectric discs of radius r = 0.37a in vacuum; the 
discs were of dielectric 

/ 1.025 -n/3/40 \ 

e = £6 -73/40 1.075 (21) 

V 1./ 

where we vary the scalar to vary the dielectric contrast. For ei, — 10, the above is the 
diagonal matrix (10, 11, 10) rotated by 30 degrees around the z-axis. 

Figure |4] shows the relative error in the frequency of the lowest several modes vs. 
the number of cells per vacuum wavelength (or c divided by the mode frequency), for 
dielectric contrast of ei, = 5, 10, 30, 100. Low contrast simulations, eb < 10, yield second- 
order error to rather high resolutions. At sufficiently high resolution, the error becomes 
first-order; this is more clearly seen in the mid-contrasts. For high contrast, ej, > 30, the 
second-order region is too small to notice, and first-order convergence is clear. 

There is a problem in examining the convergence of the fields. The fields are dis- 
continuous at the dielectric interface, and it is not obvious how best to interpolate the 
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Figure 4: (color online) For the new algoritlim presented in this paper, relative errors vs. resolution 
for mode frequencies for a 2D photonic crystal of r/a = 0.37 anisotropic discs with varying dielectric 
contrast ej. The error transitions from second-order to first-order at higher resolutions. The transition 
point occurs at coarser resolution for higher dielectric contrast. 
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fields near the interface. There is a danger, when choosing an interpolation method, 
that it might not be the best interpolation method. Therefore, we avoid the interface. 
Staying at least 3Ax (where Ax is the cell size) away from the surface, there is no serious 
ambiguity in interpolation: a simple bilinear (or, in 3D, trilinear) interpolation should 
be sufficient for errors of at least 0(Aa;^). 
We examine field-convergence in two ways. 

First, we generate a set of (evenly-distributed) points on a circle at a radius a/8 larger 
than the interface, and, at each resolution, interpolate the mode fields to those points, 
comparing against our normative values (from the algorithm of [6 ) using the ^2-norm 
over the set of points (after normalizing the entire eigenmodes). We then graph the 
relative error vs. resolution. We find (not surprisingly) that it converges at the same 
rate (ultimately 0{Ax)) as the mode frequencies (Fig. [5]). 

For £{, = 100, modes 4 and 5 have relatively large errors in the surface field; qualita- 
tively, however, the modes look more accurate than the £2 norm suggests. These modes 
place most of the field energy inside the dielectric; they resemble a pair of quadrupole 
(e.g., TE21) modes in a circular waveguide: that is, the field patterns have nearly 
(but not exactly, due to the square lattice) an azimuthal dependence cos(20 + 6q) and 
cos(20 -1-^0 + '''/2) for some angle 9o- If the dielectric were isotropic, these modes would 
be degenerate, and 9o could be chosen arbitrarily (since there is a linear combination of 
the above two terms that yields cos(26' + 6'q) for any 6q). With the anisotropic dielectric, 
the mode frequencies differ by about 0.3%, and so Oq is determined. It appears that 
the error is so high because the eigenmodes have a large error in 9q. In other words, 
the field patterns look very similar to the correct fields, except they are rotated slightly. 
This is a consequence of the difficulty of eigensolving for nearly-degenerate modes; as two 
modes approach degeneracy it becomes impossible to separate them correctly (without 
recourse to some other operator). In this light, it is not surprising that when the error 
in frequency is larger than the actual separation between the two modes, the resulting 
eigenmodes may be the wrong linear combinations of the exact eigenmodes. Indeed the 
surface error starts diminishing for modes 4 and 5 approximately when the frequency 
errors approach 0.3% (see Fig. |4|. 

Second, we generate a set of points on a circle that is a radius 3Aa; outside the inter- 
face; thus each different resolution has a different set of points; as simulation resolution 
increases, the points move closer to the actual interface. For each resolution, the fields 
are compared to our normative high-resolution simulation. In this case, the error does 
not vanish with Ax; in other words, it is 0(1) (Fig. [6|. 

This latter conclusion is disquieting for those who want particularly to measure fields 
at the interface. However, we point at that the error, although 0(1), can be quite low, 
especially for low dielectric contrast. For dielectric contrast < 15, the relative error in 
field is less than one percent. The error drops significantly as one moves away from the 
interface, and we believe it might be possible to extrapolate the fields from several cells 
away to find surface fields with vanishing error as Aa; — > 0. 

The wc07 method [SJ[7j, and indeed all the effective FDTD dielectric methods we've 
tried, have very similar error convergence, transitioning from second- to first-order at 
a resolution that decreases as the dielectric contrast increases. We believe this may 
explain a discrepancy that has puzzled us: this work and [S] see first-order error (in 
mode frequencies) for this effective dielectric, while the results of [71 show second-order 
error. The latter uses two anisotropic dielectrics, one with eigenvalues (1.45, 2.81, 4.98), 

15 



Field error for f,, = 5 

• • mode 1 : 

— mode 2 

— mode 3 - 10 

— mode 4 
+ - + mode 5 

— mode 6 " 2 1 0' 

— 0(Ax) - 

^10 
10 




10" 10" 
Cells per vacuum wavelength 



10" 



10 



10' 









• ■ mode 1 

— mode 2 

— mode 3 

— mode 4 
+ -+ mode 5 

— mode 6 

0(Ax) 

--■ 0{Ax^) 































10" 10' 
Cells per vacuum wavelength 



10' 





















• • mode 1 
' - ■ mode 2 

— mode 3 

— mode 4 
+ -+ mode 5 

— mode 6 
0(Ax) 

— - ■ 0(Ax= ) 



















10' 



10" 10" 
Cells per vacuum wavelength 



10" 



10" 
10^ 
2 10' 

m 

^10" 
10^ 
10 ' 



■ - mode 1 

— mode 2 

— mode 3 

— mode 4 
+ - + mode 5 

- - - ■ mode 6 

— 0(Ax) 
--■ O(Ax^) 












V 



























10' 



10^ 10" 
Cells per vacuum wavelength 



10" 



Figure 5: (color online) For the new method, 2D, anisotropic: relative error (in an ^2-norm) vs. resolution 
in E at points on a circle of radius a/8 outside the dielectric interface, for a 2D photonic crystal of 
anisotropic discs with varying dielectric contrast e;,. At high resolutions, the error is first-order. 
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Figure 6: (color online) For the new method, 2D, anisotropic: relative error (in an ^2-norm) vs. resolution 
in E at points on a sphere of radius (0.37 + 3Ax)a for a 2D photonic crystal of anisotropic discs with 
varying dielectric contrast ej. This error does not go to zero as Aa; — > 0. 
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Figure 7; (color online) For the new method, 2D, isotropic: relative errors vs. resolution for mode 
frequencies for a 2D photonic crystal of r/a = 0.37 isotropic discs with dielectric contrast ef, = 15 (left) 
and Ef, = 100 (right). Again, the error transitions from second-order to first-order. The transition point 
occurs at coarser resolution for higher dielectric contrast. 



and another with (8.49, 8.78, 11.52), both rotated by random orthogonal matrices. At 
contrast 11.52/1.45 « 8, we do in fact see second-order behavior up to hundreds of cells 
per wavelength, and for contrasts between other pairs (e.g., 8.49/4.98), second-order 
behavior persists up to higher resolutions than we have explored. 

In fact, practically, it must be said that the effective dielectric in [1] [7] does yield 
second-order behavior for low-contrast dielectric up to nearly the highest resolutions 
that one might practically use. Ultimately, however, it has first-order error. 

For low contrasts, the local error associated with the dielectric interface doesn't see 
to have much effect, and so the results are similar for different algorithms. 

8.2. Convergence: 2D isotropic discs 

To show that using isotropic dielectric (instead of anisotropic) does not change the 
order of error, we present frequency convergence results for the same problem in the 
previous section, except that the dielectric is replaced with an isotropic dielectric £ = £;,. 
Figure [7] shows frequency error convergence for £{, = 15 and £{, = 100. The former shows 
a gradual transition from second-order to first-order error, while the the higher contrast 
shows just first-order error. 

8.3. Convergence: 3D isotropic spheres 

Figure [8] shows frequency convergence for a 3D cubic lattice of isotropic spheres 
{r/a = 0.37) with e = 15 and e = 30. Although the computational requirements prevent 
exploration over the wide range of resolutions of the 2D simulations, the e = 15 case shows 
mostly second order behavior starting to transition to first-order, while e = 30 shows 
nearly first-order behavior. This is very similar to Fig. [4j considering the resolutions 
where they overlap. There is no reason to suspect that error convergence in 3D is any 
different from 2D. 
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Figure 8; (color online) For the new method, 3D, isotropic: relative errors vs. resolution for mode 
frequencies for a 3D cubic lattice of spheres with isotropic e = 15 (left) and e = 30 (right), using the 
algorithm of Sec. [7] For e = 15 convergence is second-order for low resolution, but gradually moves 
toward first-order; for e = 30, convergence is nearly first-order. 



S.Jf-. Comparison to wc07: 2D and 3D, aniso- and iso-tropic 

In this section we compare the "new" method recommended in this paper (Sec.[7| to 
wc07 (which uses the effective dielectric of [5]), except we improved wc07 for anisotropic 
dielectrics by using the effective dielectric of [U [7] . 

The main advantage of the new method is that it is always stable; since wc07 is usually 
stable for contrasts less than e — 30, and can be used practically, we wanted to compare 
their accuracies. The new method is generally better than wc07, by a small factor (2-10), 
for higher contrasts and resolutions. We cannot compare them for contrasts at e = 100 
because the old method becomes unstable. 

For low contrast, and low resolution, there is less reason to choose one algorithm over 
the other — for some modes one is better, for other modes the other is better. However, 
even for e = 5, the new method can yield more accurate fields, even while the frequencies 
are more or less equally accurate. Of course, there may be geometries and contrasts for 
which wc07 is better. 

Figure |9] shows the error for wc07 divided by the error for the new method, for the 2D 
square lattice of anisotropic discs (as in Sec. 8.1 ). Not only is the new method guaranteed 



stable, it performs better for medium and high dielectric contrast. 

Occasionally, for a problem of low contrast, one sees the error of wc07 plunge at some 
low resolution; the same sometimes happens for other methods as well. In this case, it 
appears that the frequency error has first- and second-order contributions: a Ax + /3Ax^. 
When a and /3 have opposite signs, there is a narrow range of Ax for which the error is 
nearly zero. However, this drop in frequency error is not reflected in the field error. An 
example of this can be seen for e = 5 in Figs. [13] and [T4j 

Figure [To] shows the error of algorithm wc07 divided by that of the new algorithm for 
surface fields on the circle at r/a — 0.37 + 1/8. In nearly all cases, the wc07 algorithm 
has higher error. The error ratio is seen to increase with the number of cells per vacuum 
wavelength. 

A similar advantage for the new method is also apparent in 3D, as shown in Fig. |11[ 
for sapphire (anisotropic, eb ~ 10) and isotropic, ei, — 15 spheres. The dielectric tensor 
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Figure 9: (color online) For wc07/new, 2D, anisotropic: the frequency error in wc07 over the error in 
the new algorithm, for modes of a 2D photonic crystal of r/a = 0.37 anisotropic discs with dielectric 
contrast ej. For e;, > 10, the new method has lower frequency error; however. Fig. |10| shows that even 
for £5 = 5, the new method has lower field error. 
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Figure 10: (color online) For wc07/new, 2D, anisotropic: the field error in wc07 over the error in the 
new algorithm, for fields on a circle at r/a = 0.37+ 1/8, a fixed distance outside the dielectric disc. The 
new method is better in almost every case, even for ej, = 5. (The simulations are the same as in Fig.^) 
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Figure 11: (color online) For wc07/new, 3D, aniso- and iso-tropic: the frequency error in wc07 over the 
error in the new algorithm, for modes of a 3D photonic crystal of r/a = 0.37 spheres of sapphire (left) 
and isotropic e = 15 (right). 
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Figure 12; (color online) For wc07/new, 2D, isotropic (inverse): the frequency error in wc07 over the 
error in the new algorithm, for modes of a 2D photonic crystal of r/a = 0.37 vacuum discs (holes) inside 
an isotropic background of e = 15. 



for the sapphire spheres was 



10.225 
-0.825 
-0.55^/3/2 



-0.825 
10.225 
0.55^/3/2 




(22) 



which has e ~ 11.6 along its c-axis, and 9.4 in the two perpendicular directions; we took 
the c-axis along y, and then rotated the dielectric by 30° about the x-axis, and then 45° 
about the z-axis. 

The advantage of the new method remains when the concavity of the dielectric inter- 



face changes: Figure 12 shows the advantage of the new method for vacuum discs in a 
background of £{,. 
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Figure 13: (color online) For wc07mod/new, 2D, anisotropic: the frequency error in wcOTmod over the 
error in the new algorithm, modes of a 2D photonic crystal oir/a = 0.37 anisotropic discs with dielectric 
contrast e;,. 



8.5. Comparison to wcOlmod: 2D anisotropic discs 

Here we compare the new algorithm to the wcOTmod algorithm (Sec.|6]), for the modes 
of a square lattice oi r /a = 0.37 anisotropic discs, as before. The new algorithm is better 
except (surprisingly, given how much worse wcOTmod is at medium contrast) for Cf, = 100. 



Figure 13 shows the ratio in frequency errors. For low contrast, we see that the error 
for wc07mod plunges within a narrow range of resolutions, due to fortuitous cancellation 
of first- and second-order error; we note that one cannot depend on this cancellation 
(it does not occur for all shapes, nor is the range of resolutions predictable). This is 
illustrated by the convergence of field, where such fortuitous cancellation is much less 
likely; there, wcOTmod has no less error. 



Figure 14 shows the ratio of the surface field error of the wc07mod algorithm to that 
of the new algorithm on a circle at r/a — 0.37 -I- 1/8. This shows that the wc07mod 
algorithm generally has higher error, and that the factor by which it is worse increases 
with the number of cells per vacuum wavelength. 

9. Summary and Discussion 

We have demonstrated a new FDTD algorithm for simulating electromagnetics in 
the presence of sharp dielectric transitions; the new algorithm is generally as accurate or 
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Figure 14; (color online) For wc07mod/new, 2D, anisotropic: the field error in wcOTmod over the error 
in the new algorithm, for fields on a circle at r/a = 0.37+ 1/8, a fixed distance outside the dielectric 
disc. 
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more accurate than previous FDTD algorithms, and it is stable at high dielectric contrast 
(unlike the algorithm of |^[7j). 

We showed how to create a stable algorithm, given the ability to form an "effective" 
(or "average") 3x3 dielectric tensor relating any three neighboring components of the 
-D-field and the ii^-field. As long as each 3x3 tensor is symmetric and positive definite 
(SPD), the algorithm will be stable. One can then try to construct SPD effective dielectric 
tensors to achieve the highest accuracy. 

We compared three different ways to compute the effective dielectric: a sometimes 
unstable method, wc07, described in [Sl[7]; a small alteration of that method to achieve 
stability, wc07mod, but still using the dielectric averaging of [1] ; and a new method based 
on symmetrizing the asymmetric effective dielectric of [6j. 

All these algorithms (except [Hj , which is unstable in the time domain) have first-order 
error in the grid-cell-size Ax; that is, the error in mode frequencies and fields at given 
points decreases ultimately as 0{Ax). However, at coarse resolutions the error decreases 
as 0{Ax'^), before transitioning to 0{Ax). For low dielectric contrast (less than about 
fO), the transition point occurs at fairly high resolution, so that for many simulations, 
these methods may be practically considered to have second-order error. 

By examining the convergence of mode fields, we showed that the fields at fixed points 
(fixed as Aa; varies) converge at the same rate as the mode frequencies. However, when 
we look at the fields at a fixed cell-distance away from a dielectric interface (so that the 
points move closer to the interface as Aa: decreases), we find that the error is 0{1) — it 
does not decrease with Ax. 

This demonstrates two important points. First, the global error is 0{Ax): the local 
error in the bulk material is O(Ax^) (due to centered differencing), is (eventually) eclipsed 
by the 0(1) local error at the interface (due to the discontinuity in fields). However, the 
interface cuts a fraction of cells scaling as 0{Ax), so the ultimate global error (e.g., in 
mode frequency) in 0{Ax). Heuristically, we can write "Error = aAx + bAx'^" where 
a comes from the interface, and b from the bulk (as well higher order terms from the 
interface). Our results show that a can be relatively small for low dielectric contrast, 
and therefore the error at low resolutions (large Aa:) appears second-order for a while. 
Sometimes a and b may be of opposite signs, in which case the error may plunge briefly 
for a small range of Aa; at the transition between second- and first-order (but only for 
frequency error; field error is not a one-dimensional quantity like frequency, and so the 
first- and second-order parts cannot fortuitously cancel). 

Although we cannot prove that the symmetric effective dielectrics used in FDTD 
simulations must yield (ultimate) first-order error, the logic of BJ strongly hints that 
should be the case. Reference ;6] achieved second-order error by demanding local first- 
order, or 0(Aa:), error at the dielectric interface — accomplished by finding an effective 
dielectric that exactly maps D to E in the limit of infinite wavelength (and planar 
interface). Unfortunately, that effective dielectric is asymmetric, and unusable in FDTD 
codes due to the instability it creates. Equally unfortunate, the symmetric effective 
dielectrics do not appear to be able to satisfy this property of mapping D to E exactly 
in the infinite wavelength limit. 

The 0(1) convergence of error at a fixed cell-distance from the surface could be a 
problem for finding surface fields with arbitrary accuracy. However, the error drops 
rapidly by several cells away from the surface, and we believe that it will be possible 
to obtain surface fields with arbitrary accuracy by extrapolating from the fields several 
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cells away from the surface. Of course, that extrapolation introduces some error, but we 
currently believe that the field error plus the extrapolation error can be made to vanish 
as Aa; — > 0. However, a much more thorough study will be needed in the future to 
determine this issue. 

It would be great if we could somehow form a symmetric effective dielectric with the 
accuracy of the asymmetric effective dielectric of [6]. Indeed, there are some degrees 
of freedom one could exploit. For example, we used a symmetric average of the S^^^ 
matrices is Sec.|4]to form S. We could have have used different non-negative coefficients 
(not 1/8) that add up to one, without affecting stability; however, the accuracy within 
uniform, anisotropic dielectric would be degraded. In principle, one can vary those 
coefficients in space; however, varying them can destroy the symmetry of S, so one 
would need some sort of global solution to attain accuracy and symmetry (not to mention 
positive definiteness) . We have made some attempts to do this, without complete success. 
There may be a way; if not, it may still be possible to achieve greater accuracy, though 
not as high as 0, while maintaining stability. 
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Appendix A. The error in the asymmetric second-order method 



To demonstrate what second-order convergence looks like, we show convergence using 
the results for the 2D anisotropic dielectric problem for €b = 10 and et, = 100 using the 
algorithm of [5] in Fig. 



A.15 



Here we compare to the Richardson-extrapolated values 
from simulations at resolutions Na = 512 and Na = 1024 (and therefore, we do not plot 
the values for Na = 1024). 



The code used to produce these graphs is open-sourced at https://github.com/ 
Ibauerca/maxwelll 



Appendix B. The effective dielectric tensor of (Kottke, 2008) is positive 
semi-definite 

If two dielectrics, and e"^, coexist in a cell in which the normal to the interface 
between dielectrics is n, then we compute the effective dielectric as follows, according to 

Rotating into a coordinate system where the first coordinate is h (the other two 
directions, perpendicular to fi, are labeled 2 and 3), we calculate and r^: 
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According to [T], r (not, e.g., e) is the quantity that should be volume-averaged. We 
perform a simple average on and to get the effective f , and then transform back 
to get the effective e. 
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For example, if the volume fractions for and are Vi and V2, respectively (V1+V2 = 1), 
then the effective dielectric is 

e^e{ViT^ +V2T''). (B.3) 

We can show that, if and are symmetric and positive semi-definite (SPSD), then 
e is also SPSD. 

We remember that a real, symmetric matrix A is positive semi-definite when, for all 
a;, Ax > 0. It follows that the sum of SPSD matrices is again SPSD. 
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Figure A. 15: (color online) For the second-order method, 2D, anisotropic: relative error vs. resolution 
for the second-order method of (which would be unstable in the time-domain), for the anisotropic 
2D photonic crystal with ei, = 10 (left) and ej, = 100 (right). Errors shown are in frequency (top), and 
in E (middle) for a circle at radius (0.37 + l/8)a, and (bottom) for a circle of radius (0.37 -I- 3Ax)a. As 
expected, global errors (top and middle) are O(Ax^), but local field errors at the interface are ultimately 
0(Ax) (as clearly seen in ej, = 100). 



28 



Whenever A is real and SPSD, there exists a matrix P such that A = P-^P{^ The 
reverse also holds: if A = P^P, then for all x, Ax = x^P'^Px = ||-Pa;|P > 0. From 
this we can conclude that Q^AQ — (PQ)^ (PQ) is SPSD for any matrix Q. 

We will prove that e = e(f) = e{ViT^ + V2T^) is SPSD, by showing that there is an 
invertible matrix F such that F^eF is SPSD. 
The matrix F 

£n2 £n3 



I 



F(£) = 



1 

'-7171 





V 



(B.4) 



is invertible, since it is diagonalizable and none of its eigenvalues (l/e„„, 1, and 1) are 
zero. The matrix product 
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is SPSD whenever e is (since x-^F-^eFx = (Fa;)^e(Fx)). Since F is invertible the converse 
also holds: if F^eF is positive semi-definite, then so is e. 

AU this means that because the e* are SPSD, F(£*)^eT(£') are SPSD. Since f = 
ViT^ + V2T^, the matrix 



Tier em 
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is the sum of two positive semi-definite matrices, which is itself positive semi-definite. 
Therefore the effective dielectric tensor, e is positive semi-definite. 

Incidentally, the matrix F plays a well-known role. If F = (Z3„, E2, E^) is the 3-tuple 
of continuous field components (at the dielectric interface), then F = F(e)E. In fact, the 
effective dielectric is derived (see [I]) by taking the perturbative formula for change in 
eigenvalue: 

d.uj^ {E\e{x) - e' {x)\E') (B.7) 

where E is the eigenmode corresponding to e{x) and E' is the (exact) eigenmode corre- 
sponding to e'{x). The above expression contains no approximation; the problem is that 
one typically doesn't know E' — usually a first order approximation is realized by approx- 
imating E' w E. In this case, the authors of 1 argued that because E is discontinuous, 
it jumps when e changes by a finite amount, and so P' ~ £' is a bad approximation (it 
has zeroth order error). It is better to approximate F' w F. Therefore, we write 

Aw^ cx {E\e{x) - e'{x)\E') = {F\T{sf[e{x) - e'(x)]F(e')|F') 



{F\r{sf[s{x)^s'{x)ne')\F). 



(B. 



^ Since A is symmetric, it has an eigendecomposition A = DV; since A is positive semi-definite, 
the eigenvalues are non-negative, and D has a real square root; therefore, A = (y/oV)'^ (y/DV). 
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To obtain approximately zero frequency shift when we substitute e' for e, we demand that 
Jy T{£)'^[e{x)—£'{x)]T{e') vanish, which leads directly to the condition that Jy T(e(a;)) = 
Ivr{e'{x)). 
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